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Abstract 

We introduce an exact reformulation of a broad class of neighborhood filters, among 
which the bilateral filters, in terms of two functional rearrangements: the decreasing and the 
relative rearrangements. 

Independently of the image spatial dimension (one-dimensional signal, image, volume of 
images, etc.), we reformulate these filters as integral operators defined in a one-dimensional 
space corresponding to the level sets measures. 

We prove the equivalence between the usual pixel-based version and the rearranged ver¬ 
sion of the filter. When restricted to the discrete setting, our reformulation of bilateral filters 
extends previous results for the so-called fast bilateral filtering. We, in addition, prove that 
the solution of the discrete setting, understood as constant-wise interpolators, converges to 
the solution of the continuous setting. 

Finally, we numerically illustrate computational aspects concerning quality approxima¬ 
tion and execution time provided by the rearranged formulation. 

Keywords: Neighborhood filters, bilateral filter, decreasing rearrangement, relative rear¬ 
rangement, denoising. 
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1 Introduction 

Let Q C R d (d > 1) be an open and bounded set, u € L°°(0) be an intensity image, and consider 
the family of filters to which we shall refer broadly as to bilateral filters, defined by 

F = cfio l n ~ u (y)) w p(\ x ~ ylMy) rf y> (i) 

where h and p are positive constants, and 

C(x) = f JC h (u(x)-u(y))w p (\x.-y\)dy 
J n 

is a normalization factor. 

Functions V/,(£) = fC(f/h) and w p are the range kernel and the spatial kernel of the filter, 
respectively, making reference to their type of interaction with the image domain. A usual choice 
for V is the Gaussian /C(£) = exp(—£ 2 ), while different choices of w p give rise to several well 
known neighborhood filters, e.g., 
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• The Neighborhood filter, see [ID], for w p = 1. 

• The Yaroslavsky filter j5T[(^5], for w p (\x — y|) = xb p (x)( y)> the characteristic function of 
a ball centered at x of radios p. 

• The SUSAN [33] or Bilateral filters [35] . for w p (s) = exp (—(s/p) 2 ). 

These filters have been introduced in the last decades as alternatives to local methods such 
as those expressed in terms of nonlinear diffusion partial differential equations (PDE), among 
which the pioneering approaches of Perona and Malik [33] . Alvarez, Lions and Morel [3] and 
Rudin, Osher and Fatemi [32] are fundamental. We refer the reader to [12] for a review and 
comparison of these methods. 


1.1 Prior work 


Neighborhood filters have been mathematically analyzed from different points of view. For 
instance, Barash [g], Elad [19] . Barash et al. [ffj, and Buades et al. [11] investigate the asymptotic 
relationship between the Yaroslavsky filter and the Perona-Malik equation. Gilboa et al. [23] 
study certain applications of nonlocal operators to hnage processing. In [[53] , Peyre establishes a 
relationship between nonlocal filtering schemes and thresholding in adapted orthogonal basis. In 
a more recent paper, Singer et al. [33] interpret the Neighborhood filter as a stochastic diffusion 
process, explaining in this way the attenuation of high frequencies in the processed images. 

From the computational point of view, until the reformulation given by Porikli [37]. their 
actual implementation was of limited use due to the high computational demand of the direct 
space-range discretization. Only window-sliding optimization, like that introduced by Weiss [33] 
to avoid redundant kernel calculations, or filter approximations, like the introduced by Paris and 
Durand [32], were of computational use. In [32]. the space and range domains are merged into 
a single domain where the bilateral filter may be expressed as a linear convolution, followed 
by two simple nonlinearities. This allowed the authors to derive simple down-sampling criteria, 
which were the key for filtering acceleration. 

However, in EH, the author introduced a new exact discrete formulation of the bilateral filter 
for spatial box kernel (Yarsolavsky filter) using the local histograms of the image, h x = /i|b p (x)> 
where B p (x) is the box of radios p centered at pixel x, arriving to the formula 


Fu(x) 



qih x (qi)ICh(u(x) - qi ), 


( 2 ) 


where the range of summation is over the quantized values of the image, qi,... ,q n , instead of 
over the pixel spatial range. In addition, a zig-zag pixel scanning technique was used so that 
the local histogram is actualized only in the borders of the spatial kernel box. 

Formula (J2[) is an exact formulation of the box filter where all the terms but the local 
histogram may be computed separately in constant time, and it is therefore referred to as a 
constant time 0(1) method. 

Unfortunately, the use of local histograms is only valid for constant-wise spatial kernels, and 
subsequent applications of the new formulation to general spatial kernels is, with the exception of 
polynomial and trigonometric polynomial kernels, only approximated. Thus, in [37] polynomial 
approximation was used to deal with the usual spatial Gaussian kernel. This idea was unproved 
in [13] by using trigonometric expansions. 
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In [IS] , Yang et al. introduced a new 0(1) method capable of handling arbitrary spatial and 
range kernels, as an extension of the ideas of Durand et al (17] ■ They use the so-called Principle 
Bilateral Filtered Image Component Jk, given by, for u(x) = q^, 

E y eAT(x) ^(g fc ~ ^(y)K(l x - yIMy) 
qk x ~ Eye^x) - u (y)H(l x - y|) 

where N(x) is some neighborhood of x. Then, the bilateral filter may be expressed as F u(x) = 
./ u( x)( x ). In practice, only a subset of the range values is considered, and the final filtered image 
is produced by linear interpolation. In this situation this filter is, thus, an approximation to 
the bilateral filter. The same authors have recently extended and optimized [60] the method of 
Paris et al. jS2j by solving cost volume aggregation problems. 

Other approaches include that of Kass et al. [26], who used a smoothed histogram to 
accelerate median filtering and proposed mode-based filters. Adams et al. [I] proposed to use 
Gaussian KD-trees for efficient high-dimensional Gaussian filtering, integrating this method with 
Paris et al. method [32] for fast bilateral filtering. They later proposed to use permutohedral 
lattice {2j for bilateral filtering, which is faster than Gaussian KD-trees for relatively lower 
dimensionality. However, both Gaussian KD-trees and permutohedral lattice are relatively less 
efficient when applied to intensity images. Ram et al. [51] used a smooth patches reordering to 
image denoising and inpainting which, when the patches are shrunk to pixels, contain similar 
ideas to the decreasing rearrangement approach. A review of some of these methods may be 
found in [26] . 


1.2 Preceding applications of functional rearrangement 

In [22], we heuristically introduced a denoising algorithm based in the Neighborhood filter for 
denoising images produced as time-frequency representations of one-dimensional signals, for 
which the large computational effort made useless other denoising methods based on PDE’s. 
The main observation was that the Neighborhood filter can be computed using only the level 
sets of the image, instead of the pixel space, producing a fast denoising algorithm. 

In [23] . this idea was rendered to a rigorous mathematical formulation through the notion 
of decreasing rearrangement of a function. In short, the decreasing rearrangement of u, denoted 
by «*, is defined as the inverse of the distribution function m u (q) = |{x E f2 : u(x) > g}|, see 
Section 2 for the precise definition and some of its properties. 

Realizing that the structure of level sets of u is invariant through the Neighborhood filter 
operation as well as through the decreasing rearrangement of u. allowed us to rewrite (PQ), for the 
homogeneous kernel w p = 1 , in terms of the one-dimensional integral expression 


NF* u(x) 


JCh(u(x) - u*(s))u*(s)e/s 
/J n| JC h (u(x) - u t (s))ds 


(3) 


winch is computed jointly for all the pixels in each level set {x : u(x) = q}. 

Perhaps, the most important consequence of using the rearrangement was, apart from the 
large dimensional reduction, the reinterpretation of the Neighborhood filter as a local algorithm. 
Indeed, notice that since u * is decreasing, we have 


u*(f) — u*(s)| < h- if \t — s\ < c(h), 
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for some continuous function c. Thus, the range kernel of the pixel-based filter (U), whose support 
in the pixel space may be disconnected, is transformed into a local kernel in its rearranged version 


Thanks to this we proved, among others, the following properties for the most usual nonhnear 
iterative variant of the Neighborhood filter: 

• The asymptotic behavior of NF* as a shock filter of the type introduced by Alvarez et al. 
(1], combined with a contrast loss effect. 

• The contrast change character of the Neighborhood filter, i.e. the existence of a continuous 
and increasing function g : M —> K such that 

NF u(x) = NF» u(x) = g{u{x)). 


1.3 Objective of the article 


In this article we extend the use of functional rearranging techniques, as introduced in [23], to 
bilateral filters which, in the discrete case and for the spatial box window, coincides with the 
bilateral filter reformulation given by Porikli m- 

We provide a rigorous mathematical ground which allows to interpret the discrete formulas 
of [37], and their extension to general spatial and range kernels, as the discrete counterpart of 
continuous pixel-space formulations. 

The formulas we obtain for general bilateral filters allows to reinterpret these filters from 
the range space point of view, as a natural extension to Porikli’s discrete model. In addition, 
these formulas may be computationally competitive when a high degree of approximation to the 
pixel-based formulation is required. 

This work may be considered as a non trivial extension of the rearranging techniques in¬ 
troduced in [23] for the shnpler Neighborhood filter. Thus, even if the spatial kernel w p is 
non-homogeneous, we may still use our approach by introducing the relative rearrangement of 
the spatial kernel with respect to the image, see Section 2 for definitions. 

Using this tool, we may express the general bilateral filter jT]) in terms of the one-dimensional 
integral expression 


F*u(x) 


JT £h(u{x) - u*(s))w p { |x - 
/J n| ICh(u(x) - u*(s))w p (|x 


j )*u(s)u*(s)ds 
- j)*u(s)ds 


(4) 


where v* u denotes the relative rearrangement of v with respect to u. 

To gain some insight into formula Q, let us consider a constant-wise interpolation of a 
given image, u, quantized in n levels labeled by q ,, with max('u) = q\ > ... > q n = 0. That is 
u(x) = £r=i q z XE, ( x ), where Ei are the level sets of u, 


Ei = {x G ft : u(x) = qi}, i = l,...,n. 


Similarly, let w p be a constant-wise interpolation of the spatial kernel quantized in m levels, rj, 
with max(wp) = ri > ... > r m = min(t/7 p ) > 0. For each xgO, consider the partition of E t 
given by Fj(x) - {y £ £, : «; p (|x — y|) = rj}. Then, we show hi Theorems [T] and |3] that for each 
x e E k , k = 1,... ,n, 


F u(x) = F, u(x) 


E”=i £h(qk - qi)Wim(x)gi 

Er=l IC h (q k -qi)Wim(x) ’ 


(5) 
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where W ITn (x) = ( x )|» with | ■ | denoting set measure (number of pixels). That is, 

|F) (x)| is the number of pixels in the gj-level set of u that belongs to the r^-level set of w p . 

Observe that, like for the Neighborhood filter, the range kernel is transformed into a local 
kernel, independent of the pixel location, and it is therefore computed once and for all, explaining 
the large gain in computational cost, as already observed in [37] . 

In fact, notice that in the particular case in winch w p is taken as a box window we have that 
Wim(x) is just the local histogram, h x , used in Porikli’s formula, so ((21) and (HJ) coincide in this 
case. However, the rearranged formula (2J) is valid for any other shape of the spatial kernel. 

The results in this article may be summarized as follows: 

1. The rearranged formulation © is exact and has no restrictions on the spatial or range 
kernels shape (Theorem [I]). In fact, our results are not limited to the special functional form of 
the spatial kernel, w p = u> p (|x —y|), that we considered along this article. One can easily modify 
our arguments to deal with more general kernels of the type w p = w p (x. y), thus including other 
well known filters like the median filter or the cross-joint bilateral filter TH] 135] • We preferred 
to limit ourselves to bilateral type filters for the sake of clarity, and leave the more general case 
to future developments of the work. 

2. When restricted to a spatial box kernel, the relative rearrangement term of formula flU) is 
just the local histogram, like in Porikli’s discrete formulation. For more general spatial kernels, 
the relative rearrangement may be interpreted as an averaged local histogram evaluated on the 
kernel level sets, extending in this way the notion of local histogram and the discrete formulation 
of [57] . 

3. We establish the relationship between the discrete and the continuous image setting of the 
rearranged version of the bilateral filter, proving the convergence of constantwise interpolators 
to their continuous counterpart (Theorems [21 (31) ■ 


1.4 Outline 

The plan of the article is the following. In Section 2 we recall the notions of decreasing and 
relative rearrangements and motivate the deduction of the rearranged formula Q under some 
restrictive regularity assumptions. 

In Section 3, we remove those assumptions and establish the equivalence between the usual 
pixel-based expression of the filter (PQ) and its rearranged formulation dU). hi addition, we provide 
a fully discrete algorithm to approximate by constant-wise functions the filter F*u(x) given by 
@, and thus the original equivalent filter Fu(x) given by (JTJ). We also prove the convergence 
of this discretization to the solution of the continuous setting. In Section 4, we give the proofs 
of these results. 

In Section 5 we illustrate the performance of the rearranged filter formulation in comparison 
to the brute force pixel-based implementation and to other state of the art filters: the 0(1) 
method of Yang et al. [531, and the permutohedral lattice method of Adams et al. [T]. 

Finally, in Section 6 we give the Summary. 
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2 Neighborhood filters in terms of functional rearrangements 

2.1 The decreasing rearrangement 

Let us denote by |£'| the Lebesgue measure of any measurable set E. For a Lebesgue measurable 
function u : 0 —> R, the function rn u (q) = |{x £ 0 : u(x) > <j}| is called the distribution 

function correspondmg to u. 

Function m u is non-increasing and therefore admits a unique generalized inverse, called the 
decreasing rearrangement. This inverse takes the usual pointwise meaning when the function u 
has no flat regions, i.e. when |{x £ ft : u(x) = q}\ =0 for any q £ R. In general, the decreasing 
rearrangement u t : [0, |0j] —> R is given by: 

{ ess sup{w(x) : x £ 0} if s = 0, 

inf {q £ R : m u (q ) < s} if s £ (0, \Q\), 

essinf{u(x) : x £ 0} if s = |fl|. 


We shall also use the notation 0* = (0, |0|). Notice that since «* is non-increasing in 0*, it is 
continuous but at most a countable subset of 0*. In particular, it is then right-continuous for 


all t £ (0, |f2|]. 


The notion of rearrangement of a function is classical and was introduced by Hardy, Little- 


wood and Polya [25] . Applications include the study of isoperimetric and variational inequalities 
[55] [7, [3D], comparison of solutions of partial differential equations [35] [5] [37] [13] [5], and others. 

We refer the reader to the textbook [27] for the basic definitions. 

Two of the most remarkable properties of the decreasing rearrangement are the equi-measurability 
property 


r rl^l 

/ f(u(y))dy = / f(u*(s))ds, 

Jq. J o 

for any Borel function / : R —» R+, and the contractivity 


for u,v £ L P (S7), p £ [l,oo]. 


\u. 


~ v *I|lp( 0,) < ll M - w llLP(fi)i 


( 6 ) 


2.2 Motivation 

Apart from the pure mathematical interest, the reformulation of neighborhood filters in terms 
of functional rearrangements is useful for computational purposes, specially when the spatial 
kernel, w, is homogeneous, that is w p = 1. In this case, it may be proven [23] that the level sets 
of u are invariant through the filter, i.e. u(x) = u(y) implies F(u)(x) = F(u)( y), and thus, it 
is sufficient to compute the filter only for each (quantized) level set, instead of for each pixel, 
meaning a huge gain of computational effort. 

For non-homogeneous kernels the advantages of the filter rearranged version are kernel- 
dependent, and in any case, the gain is never comparable to that of homogeneous kernels. The 
main reason is that the non-homogeneity of the spatial kernel breaks, in general, the invariance 
of level sets through the filter application. 

In the following lines, we provide a heuristic derivation of the Bilateral filter rearranged 
version, as first noticed in [23]. We shall prove hi Theorem [U that the resulting formula is vahd 
in a very general setting. 
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Under suitable regularity assumptions, the coarea formula states 

[ 0 (y)|Vu(y)|dy = f f g(y)dT(y)dt. 

J Q, J —oo J u=t 


Taking g(y) = IC h (u(x.) - u{y))w p (\x - y|)u(y)/|Vu(y)|, and using u(x) € [0, Q] for all x e ft 
we get, for the numerator of the filter dTJ) , 

J(x) : = J £h(u(x) - u(y))w p (\x - y\)u{y)dy 

Assuming that m* is strictly decreasing and introducing the change of variable t = u t (s) we find 

du#(s) 


f l“l , 

/(x) = - / ACh(u(x) - u*(s))m*(s)- 
Jo 


ds 


/ 

J U=1 


w p(\* ~ y|) 


dT(y)ds 


=u.(s) |Vu(y)| 

rW 

= / /Cfc(w(x) - u*(s))u*(s)w p (\x - -\)* u {s)ds. 
Jo 


(7) 


Here, the notation v tu stands for the relative rearrangement of v with respect to u which, under 
regularity conditions, may be expressed as 


v* u (s) = 


L 

L 


yjy) 

Ju=u,( S ) |VM(y)l 


u=u,(s) 


|Vw(y) 


dT(y) 

■dT(y) 


( 8 ) 


see the next subsection for details. Transforming the denominator of the filter, C(x), in a similar 
way allows us to deduce 


Jo ^ 1 £h(tt(x) - u*(s))w p (|x - -| )* u (s)u*(s)ds 


Fu(x) = 


,/f £ h (u(x) - u*(s))w p {|x - -\) tu (s)ds 


(9) 


2.3 The relative rearrangement 


The relative rearrangement was introduced by Mossino and Temam [2b] as the directional deriva¬ 
tive of the decreasing rearrangement. Thus, if for u € L l {Q.) and v 6 L p (fl), with p £ [1, oo], we 
consider the function 

<p{s)= / w(x)dx+ / (w|«=u.( a ))*(^)d<T, 

Ju>u,(s) JO 

then the relative rearrangement of v with respect to u, v* u , is defined as 


v* u := - 7 -<p £ 
as 
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This identity may also be understood as the weak L P (G*) directional derivative (weak* L°°(S2*), 
if p = oo), 

(u + tv)„-u* 

v* u = Inn---. (10) 

Under the additional assumptions u G VU 1,1 (f2) and |{y G 0 : Vu(y) = 0}| = 0, i.e. the non¬ 
existence of flat regions of u, the identity (JBJ is well defined and coincides with (fTTTf) . In this case, 
the relative rearrangement represents an averaging procedure of the values of v on the level sets 
of u labeled by the super level sets measures, s. 

When formula (J8J) does not apply, that is, in flat regions of u, we may resort to identity (flTTj) 
to interpret the relative rearrangement as the decreasing rearrangement of v restricted to such 
flat regions of u. 

After the seminal work of Mossino and Temam [2(5], the relative rearrangement was further 
studied by Mossino and Rakotoson |31j and applied to several types of problems, among which 
those related to variable exponent spaces and functional properties [2Ql EH SS 3D|, or nonlocal 
formulations of plasma physics problems [151 (TB]. In the rest of the article, we shall make an 
extensive use of the results given in the monograph on the relative rearrangement by Rakotoson 

m- 


2.4 Example: Rearrangement of constant-wise functions 

Consider the constant-wise functions u.v : [0,13] —» M given in Fig. [T] (a). Writing max(u) = 
5 = gi>...><3 , 6 = 0 = min(u), we may express u as u(x ) = Yl*i=i QiXEi( x ), where Ei are 
the level sets of u, E\ = (10,11], Eo = (8,10], etc. Then, the decreasing rearrangement of u is 
constant-wise too, and given by 

n 

u*(s ) = '^2q iX i i (s), 

i =1 

with Ii = [ai_i, di) for i = 1,... ,6, and with ao = 0, ai = \Ei\ = 1, ai = |Ei| +1£2| = 3,.. .,06 = 
Y^=i \E>\ = 1^1 = 13- The corresponding plot is shown in Fig. Q] (b). 

In Fig. [T] (c) we show the graphs {E t , v(Ei)}f =1 transported as it was done in the step before 
to construct u«. For instance, the highest level set of u, E\ = (10,11], was transported to [0,1]; 
E 2 = (8,10] to (1,3], etc. Thus {E\,v(Ei)} = {(10,11], {2}} is transported to {[0,1], {2}}; 
{E 2 ,v(E 2 )} = {(8,10], {3}} is transported to {[1,3], {3}}, etc. 

Finally, to obtain the decreasing rearrangement of v with respect to u, u* u , we rearrange 
decreasingly the restriction of v to the level sets of u, E t , as shown in Fig. [T] (d). 


3 Main results 

In this section we provide analytical results showing that, given an image, u, a spatial kernel, 
w p , and a range kernel, /Q ( , all bounded in L°°, we may approximate the filtered image v = F u 
by a sequence of constant-wise functions, v n , which are obtained by filtering a constant-wise 
approximation of u, u n , through a discrete filter, F m , involving a constant-wise approximation 
of w p , w pTn . 

In addition, and fundamentally, we show that the filter and its approximations may be 
obtained using the rearranged version ([D]). and give explicit formulas for their computation. 
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(a) Functions u (blue) and v (red) (b) Decreasing rearrangement of u, 

u. 



(c) u, (blue) and transported v (d) The relative rearrangement v, u 
(red) 

Figure 1: Example of construction of the relative rearrangement, (c) shows u* and v transported as by 
the displacement of the level sets of u in the construction of u „. (d) shows the decreasing rearrangement 
of v restricted to the level sets of u, that is v* u . 


The main assumption we implicitly made for the heuristic deduction of formula ({5J) is the 
condition |{y G ft : Vu(y) = 0} | =0, i.e. the non-existence of flat regions of u, which gives sense 
on one hand to formula ([ 8 ]) . and on the other hand, allows us to obtain the strictly decreasing 
behavior of u*, which justifies the change of variable in (171). 

Our first result is that formulas (UJ) and © for the pixel-based filter and its rearranged 
version are equivalent under weaker hypothesis. Due to the nature of our application, we keep 
the assumption on the boundedness of u and w in L°°, although these can be also weakened. 
Theorems |T] and [2] are proven in Section g] We use the notation M + = {t G R : t > 0}. 


Theorem 1 Let Ll C R d be an open and bounded set, d > 1, K, 6 L°°(R,R + ) and w p G 
L°° (IRLf,!Lf). Assume that u G L°°(D) is, without loss of generality, non-negative. Consider 
F n(x) given by (Q3) and 


F*u(x) 


/o ^ 1 ICh{u(x) - u*(s))u> p (|x - -|)* u (s)u*(s)ds 


/ c | n ' £h(u(x) - u*(s))w p { |x - -| )* u (s)ds 


( 11 ) 
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Then F* u(x ) = F u(x) for a.e. x £ ft. 

The next step towards the definition of a fully discrete algorithm to approximate the filter 
F* given by (fTTj). and thus the original equivalent filter F given by (U), is showing that we may 
approximate u and w p by constant-wise functions u n and w pTn such that, as m , n —> oo, 

F™« n (x) ->■ F* u(x) 


where F™ is the discrete version of F*, 

pm u ( x ) = Jo 01 **M X ) - ^(s))w p7n (|x - -\) w (s)u*(s)ds ^ 

if ICh{u(x) - u*(s))w pTn (\x - -| )* u (s)ds 

Theorem 2 Assume the conditions of Theorem^ and suppose that u has a finite number of flat 
regions. Then, there exist sequences of constant-wise functions u n , w pTn , with a finite number 
of flat regions, such that u n —> u strongly in L°°(fl), w pTn —> w p strongly in and 

Ff L u n —> F#u a.e. in ft and strongly in L°°(ft). (13) 

Remark 1 Theorem [2j may be extended to the case in which u has a countable number of flat 
regions. However, in this case the approximation sequence of constant-wise functions u n has also 
a countable number of levels. Since our aim is providing a finite discretization for numerical 
implementation, such a sequence is not appropriate. 

Theorem [2] gives utility to Theorem [3} in which we produce the fully discrete formula actually 
used for the numerical implementation. 

The main difficulty in computing formula (|l2|) is the determination of the relative rear¬ 
rangement. However, in the case of constant-wise functions with a finite number of levels this 
computation is simplified thanks to identity (11 ()f) . which may be easily applied to this situation, 
as shown in [2E1 Th.7.3.4]. 

In few words, for the case of constant-wise functions u and v, the relative rearrangement v tu 
may be computed as the decreasing rearrangement of v restricted to the level sets of u. 

Theorem 3 Let u G L°°(fi) be a constant-wise function quantized in n levels labeled by q{, with 
max(u) = q\ > ... > q n = 0. That is u(x) = 9iX£i( x )> where Ei are the level sets of u, 

Ei = {x e fl : u(x) = qi}, i = 1 ,..., n. 

Similarly, let w p £ L°°(R+) be constant-wise and quantized in m levels, r ? , with ma x(w p ) = 
ri>...> r m = min^p) > 0. For each x £ ft, consider the partition of Ei given by Fj(x) = 
{y £ Ei : ie p (|x — y|) = rj}. Then, for each x £ Ek, k — 1,... ,n 

F u(x) = ^i=l K-h(qk — 

* 1 ; zr = l IC h (q k -qi)W im (x) 



with W im (x) = E^iG'l F j( x )l- 
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3.1 Examples 


As it clear- from formula (fH|) , the main difficulty for its computation is the determination of the 
measures of the sets Fj (x), which must be computed for each x 6 fi. 

The formula, provides the algorithm complexity. If N is the number of pixels, then the 
complexity (without any further code optimization) is of the order 0(Nnm ), where n is the 
number of levels of the image and m is the number of levels of the spatial kernel. Let us 
examine some examples. 

The Neighborhood filter. In this case, w p = 1, and therefore j = 1 and Fj(x) = F, is independent 
of x for all * = 1,... ,n. Thus, formula (fPIT) is computed only on the level sets of u, that is, for 
all x G Ek 




YJi=\K h {qk - qi)\E l \q l 

U=l^h(qk-qm\ 


In this case, the complexity is of order 0(n 2 ). 

The Yaroslavsky filter. In this case, w p (|x — y|) = Xb p (x) (y)i and therefore there are only two 
levels n = 1, r '2 = 0 of w p corresponding to the sets 


-P’i(x) = {y e Ei : |x - y)| < p}, 
F 2 ‘(x) = {y 6 Ei : |x - y)| > p}, 


Thus, formula ([H|l reduces to: for each x 6 Ek, k = 1,... ,n 

p* / \ = Etl^hjqk - qj)\F{(x)\qi 

h [ ’ £?=!«»-®)|Ff(x)| 

The complexity is then of order O(Nn). 

The Bilateral filter. In this case, w p (|x — y|) = exp( —|x — y | 2 /p 2 ) and therefore there is a 
continuous range of levels of w p . However, for computational purposes the range of w p is 
quantized to some finite number of levels, determined by the size of p. Thus, the full formula 
(HI must be used in this case. The resulting complexity is of order O(Nnm). 

The corresponding complexities for the pixel based formulation is, at least, of the order 
0(Nn/i), where p. denote the discrete size of the spatial kernel. However, meanwhile for the 
rearranged version the range kernel is computed only once, for the pixel-based version it must 
be computed for each pixel. 

Observe that more optimized code and a further reduction of algorithm complexity may be 
reached by standard zig-zag techniques of window scanning, see [37]. 


4 Proofs 

Proof of Theorem |T[ Let / € L°°(E) and b £ L°°(0). We start showing 

/ f(u(y))b(y)dy = / /(u*(s)) 6 * u (s)ds. 

./fi Jo 

Consider the sets of flat regions of u and rt*, 

p = IJ Pu p i = {y e n: «(y) = &}, 

ieD 


(15) 


(16) 
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and P* = U i££)P*, with , P* = {s G : u t (s) = <jj}, where the subindices set D is, at most, 
countable. According to [5B1 Lemma 2.5.2], we have 


rW r 

/ f (u*(s))b tu (s) = / /(u*(m u (u(y))))&(y)dy 

Jo Jn\p 

+ Y I M bi( h i)(y) b (y)dy, 

where bi = b\p { , /q(s ) = /(u*(s- + s)) for s G [s', s'-) := P*, and 

hi(m b .(bi( y))) if y G P,\Q\ 


M bi (hi)( y) = 


1 

\q~\ j , hi(s)ds, if y G Q), 


(17) 

(18) 


where Q l = U j^D'Q), with Q* the flat regions of b t and [o'. a") := Q 3 *. In (fT7]) . since the 
functions u* and m u are strictly decreasing and inverse of each other in the set S ~l\P, we obtain 


[ f(u*(mv(u(y))))b(y)dy = f f (u(y))b(y)dy. 

Jn\p Jn\p 


'n \p Jq\p 

In the flat regions of u and u* we have, on one hand, 


f f(u{y))b(y)dy = V / f{u(y))b{y)dy 
J P i€D JPi 

= Y /(®) / b (y) d y- 


And, on the other hand, since fi t -(s) = /(u*(s' + s)) = f(qi) for s G P*, 


Y I M bi(hi)(y)Ky)dy 

i( zD' ,Pi 


= Y\ ( h i( m bMy))) b (y)dy 

\ Jp W 


+ v 

j&D' 



b (y)dy 




b (y)dy. 


Therefore, both sides of (fT5j) are equal. 

Finally, for fixed x G set b(y) = w p (x, y) and first, /(t) = ICh(u(x) — t)t, for t > 0 to 
obtain, using the identity (fTH]) . the equality between the numerators of (HI) and (fTIjl. and second, 
/(f) = Kh(u(x ) — t) to obtain the equality between the denominators of those expressions. □ 
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Remark 2 As deduced in the proof, identity lilty follows from Lemma 2.5.2]. In fact, a 
little more than \15\) may be obtained. Let f,b G L°°(f2). Then, if f is constant in the flat 
regions of u, that is f\p i = fi = const, then 


r i“i r 

/ f(s)b* u {s)ds = / f(m u {u(y)))b{y)dy 
./ o Jn 

+ ]T/(ft) I b(y)dy. 


(19) 


Proof of Theorem [3} We split the proof in several steps. 

Step 1. We use the construction of the sequence of constant-wise functions u n given in jSS], Th. 
7.2.1]. In our case, the construction is simpler because u has a finite number of flat regions, 
implying that u n has a finite number of levels. 

In any case, this construction is such that u n —> u a.e. in Q and strongly in L°°(Q), 
and w, Un —> w* u weakly* in L°°(D*). Besides, due to the strong continuity of the decreasing 
rearrangement we also have («„)* —> it* strongly in L°°(fl*). Therefore, we readily see first 
that 

F£u n (x) —> F£u(x) for a.e. x G f2, 
and then, due to the dominated convergence theorem, 


F)*u n —> F£u strongly in L°°(Ll). 


( 20 ) 


Step 2. We consider a sequence of constant-wise functions with a finite number of levels, w m , 
such that Wm — > w strongly in L°°(S2 x f l). Due to the contractivity property of the relative 
rearrangement, see [29] , we also have, for a.e. x G O, w m {x, •)*„ —» w p (x, •)*„ strongly in L°°(D), 
for any v G L°°(D). Thus, as m —> oo, 


F^ m u n (x) = F£u n (x) for a.e. x G ft, 
and, again, the dominated convergence theorem implies 


F h,m u n -y F h u n strongly in L°°(ft). ( 21 ) 

Step 3. In view of (f2U|) and (j2T|). we have 

I F h,m u n — F h U I —\Fh,m U n ~ F h U n\ 

+ |Ffr Un — Fj£u\ —t 0 

as m —^ oo and n —> oo, so (fETI) follows. □ 

Proof of Theorem © Since u is constant-wise, the decreasing rearrangement of u is constant- 
wise too, and given by 

n 

u*(s) = 5^ftX/i(s), 

i=l 

with Ii = a*) for i = 1,... ,n, and a 0 = 0, ai = |£i|, a 2 = |£i| + \E 2 \,. • •,«« = DILi \ E i\ = 

|ft|. It is convenient to introduce here the cumulative sum of sets measures 

i 

cum(£' o ,0) = 0, and cum (E 0 ,i) = ^ \Ek\, 

k =l 
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for i = 1, where the symbol o denotes the summation variable. Thus, a t = cum(£ 0i ?'). 

We consider, for fixed x G and t > 0, the function H (y) := u( y) + tw p (|x — y|). Since 
both u and w p are constant-wise with a finite number of levels, we have, for t small enough 

Qi+l < Qi T t r j ^ Qi—h 

implying that each level set of H is included in one and only one level set of u. Thus, we have 
the orderings 

• For each j, if i\ > i 2 and y G F* 1 (x), y G Fj 2 (x) then H{ y) < y). 

• For each i, if ji > j 2 and y G F^(x), y G F] 2 (x) then H( y) < H{ y). 

With these observations we may compute the decreasing rearrangement of H as follows. For 
instance, for sG/i = [0, |£i |) we have (omitting x from the notation F/ (x)) 

gx+tr! if s G [0, \Fl\), 

Hm{s) = J 9i+*r 2 if s G [|i r 1 1 |, \Fi \ + 1^2 |), 


Notice that 


qi + tr m if s G [cum(F 0 1 ,m - 1), cmn(F 0 1 , m)). 

[0,|F?|) = [cum(F’ o 1 ,0),cum(F o 1 ,l)), 

[|*il» I*! 1 1 + l* 2 l) = [cum(F 0 1 , l),cum (-F 0 \2)), 

[ciunfF,, 1 , m — l),cun i(F^,m)) = [cmn(F 0 1 ,m - 1), |£i|), 

, in general, we may write for i = 1 ,..., n and j = 1 , ,m 


where Jj(x) := 
Finally, since J‘ 


H*{s) =qi + trj if s G Jj(x), 

6 *_ 1 (x), 6 *(x)^, with /j*(x) := aj_ 1 + cum(F*(x), j). Observe that b l m (x) = a*, 
x) C Ei we have for s G Jj(x) 

#*(«) - u*(s) . 1 . ,, j\ / \ 

- j -= rj, implymg w p (|x - -|)* u (s) = rj. 


t 

We are now fix disposition to compute formula (fTT|l . For x G Ek, 

r\Q\ 

/ IC h (u(x) - u*(s))u*(s)w p (\x - •|)* u (s)ds 
Jo 

n 771 


II III 

= M?* - qi)qirj\Jj{x) 

i =1 j =1 


In a similar way we obtain 

r\n\ 

C(x) = / IC h (u(x) - u*(s))w p ( |x - • \)* u (s)ds 
Jo 

71 771 

= IZ S Kh '( qk - q ^ r J 1 (*) 1 

i=i j=i 

and therefore, using the definition of the sets Jj(x) we deduce (TTT[) . □ 
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5 Experiments 

We conducted several experiments on standard natural images to check the performance of the 
discrete formula ([Till in comparison to well known state of the art denoising algorithms: Yang’s 
et al. Real Time 0(1) Bilateral Filtei@ [35], and the Permutohedral Lattice Fast Filtering 
[2- The aim of our experiments was to investigate the quahty of the approximation of these 
algorithms to 

1 . the ground truth, and 

2 . the exact pixel-based bilateral filter. 

We also compared execution times as delivered straightly from the available codes. Notice 
that execution time depends on code opthnization and therefore a rigorous study of this aspect 
requires some kind of code normalization which was out of the scope of our study. 

In both experiments we used four intensity images of different sizes to compare execution 
times and peak signal to noise ratio (PSNR): Clock (256 x 256), Boat (512 x 512), Airport 
(1024 x 1024), and Still life (2144 x 1424). The first three images were taken from the data 
base of the Signal and Image Processing Institute, University of Southern California, while the 
fourth was courtesy of A. Lubroth. 

We apphed the denoising filters to the test hnages corrupted with an additive Gaussian white 
noise of SNR = 10, according to the noise measure SNR = a(u)/a(u), where a is the empirical 
standard deviation, u is the original image, and u is the noise. We repeated the experiments 
described in this section for higher levels of noise. Since the results were qualitatively similar, 
we omitted them for the sake of brevity. 

We considered different spatial window sizes determined by p, with p = 4 , 8 , 16, 32. For 
the Yaroslavsky filter, the spatial kernel is the characteristic function of a box of radius p, while 
for the bilateral filter a Gaussian function is taken in the experhnents. 

The range size of the filter was taken, throughout all the experiments, as h = p which, 
according to [ 11 ] , is the regime in which the corresponding iterative filter behaves asymptotically 
as a Perona-Malik type filter. The shape of the range filter is always a Gaussian. 

The discretization of the pixel-based and the rearranged version of both filters was imple¬ 
mented in non-optimized C++ codes in which the main differences are those corresponding to 
formulas (]?3]) and (fMI) . Execution time was measured by means of function clock. All the 
experiments were run on a standard laptop (Intel Core i7-2.80 GHz processor, 8 GB RAM). 

5.1 Yaroslavsky filter 

The performance of formula (114[) for filter computation strongly depends on the functional form 
of the spatial kernel w p . In the case in which w p is the characteristic function of a box of radius 
p, that is 

F = 77TT / M“( x ) - u(y))u(y)dy, (22) 

C(x) J Bp (x) 

1 C++ code downloaded from 
http: / / wwu. cs. cityu. edu. hk/~qiyang 

^C++ code downloaded from 
http://graphics.Stanford.edu/papers/permutohedral 


15 





G. Galiano and J. Velasco 


Algorithm 1 Yaroslavsky filter (rearranged version) 

1 : * Algorithm according to formula ()24|) , with = k* 

2: for each pair (i,k) of image quantization levels do 
3: ICh(i,k ) = exp(-((k -i)/h) 2 ) 

4: (* Main loop *) 

5: for each pixel, x, of the image do 

6: Wn G- 0 

7: if x is the first pixel then 

8: for each y G B p (x) do 

9: if u(y) = i then 

10: Wii(x)+ = 1 

11: else 

12: for each y G B p (x) do 

13: actualize IFa(x) adding/deleting new/old bins w.r.t. the previous box 

14: for each quantization level, i do 

15: numerator + = ICh(i,u(x)) * Vkii(x) * i 

16: denominator + = ICh(i,u(x)) * H / ii(x) 

17: uFiltered(x) = numerator / denominator 


the number of level sets of the spatial kernel is just two, and the rearranged version requires 
only the computation of the measure of the sets F£(x) for each pixel x, that is, counting among 
all pixels in the i —level set of u how many of them belong to the spatial box. 

In the discretization of the pixel-based filter, for each pixel x we compute u in the square 
neighborhood, iV(x), centered at x of radius 2/i, containing (4/i + l) 2 pixels (close to the border, 
the image is extended by zero). Then, we approximate (j%2|) as 


E y eiv(x) ^(“W - “(y)My) 
E ye *(x)M«(*)-«( y)) 


(23) 


For the rearranged version we compute, for each pixel x G Ek, k = 1,..., n, the number of pixels 
of the spatial box which belong to the different level sets of u, Wn (x) = (M / n(x),... , VF n i(x)). 
That is, VF;i(x) is the number of pixels of E l which are present in the box N(x). Then, we take 


F* u(x) & 


IXi £h(<ik ~ gi)^ i (x)g i 

Ei=i £h(«*-«)Wii(x) ’ 


(24) 


where q = (qi,... ,q n ) are the quantization levels of u. Observe that since the reaiTanged version 
transforms the range kernel /C^(u(x) — u( y)) into, tCh(qk ~ q), this term may be computed and 
stored outside the main loop running over all the pixels, see Algorithm [T] 

Thus, in (HU), only M r (x) and qk must be actualized for each pixel, establishing an important 
difference with respect to the pixel-based version (f23|) . in winch both the spatial and range kernels 
are actualized. In addition, following EH. we perform the actualization of the measures M r (x) 
using previously computed information of neighbor boxes (zig-zag scanning technique). 
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Algorithm 2 Bilateral filter (rearranged version) 

1 : * * Algorithm according to formula (]26|) , with = k* 

2: for each pair (i, k) of image quantization levels do 
3: ICh{i,k) = exp(-((k - i)/h) 2 ) 

4: for each j of spatial kernel levels do 
5: r(j) = exp(-(j/p) 2 ) 

6: for each pixel y in a generic Gaussian window do 
7: if exp(-(y/p) 2 ) « r (j) then 

8: PixelLevel(y) j 

9: (* Main loop *) 

10 : for each pixel, x, of the image do 
11: countBins <— 0 

12 : for each pixel, y, of Gauss, window centered at x do 

13: count B in s(u(y), PixelLevel(x. — y))+ = 1 

14: for each quantization level, i do 

15: W; m (x) f- 0 

16: for each spatial kernel level, j do 

17: Wi m (x)+ = r(j) * count B ins (i, j) 

18: numerator + = fCh(i , u ( x )) * Vkj m (x) * i 

19: denominator + = ICh(i, u(x)) * Wj m (x) 

20: uFiltered(x.) = numerator / denominator 


5.2 Bilateral filter 

If the spatial kernel w p has a more complex profile than the characteristic function, e.g. the 
Gaussian function, then the rearranged version requires the computation of, for each pixel x, 
the measure of the sets Fj (x) for all the level setss, 1 < j < m, of some discretization of w p . For 
instance, for p = 4,8, 16, 32, the corresponding number of the Gaussian function quantization 
levels arising from double precision arithmetic are m = m(p), with m(4) = 42, m(8) = 135, 
m(16) = 457, and m(32) = 1621. 

The approximations to the pixel based and the rearranged version are, respectively, 

X^yeN(x) ^Wi('u(x) — rt(y))rcp(|x — y|)«(y) 

E yeiV ( x )^(«(x)-u(y)K(|x-y|) ’ 1 ° } 

and 

F 7/1-vl ~ E;=l fchjqk ~ 

* U ~ u=l Kh{qk-qi)Wim{x) ’ 

with Wi m (x.) = 1 r j |F'(x)|, where, for x 6 Ffc, |F)(x)| is the number of pixels of F, which 

belong to the j— level set of w p . 

The bilateral filter may be accelerated by manipulating the quantization levels of the image, 
and of the spatial range. As shown in jH] for the Yaroslavsky filter, the reduction of the image 
quantization levels leads to poor denoising results. However, we checked that a similar restriction 


( 26 ) 
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applied to the spatial kernel reduces the execution time while conserving good denoising quality. 
In our experiments we tested this ansatz with a fixed number of levels, m = 20, for all the 
images and p— values. 

5.3 Other algorithms and discretization details 
The algorithm introduced by Yang et al. [32] is the following. Let 

{qi,---,qa} C {?i,..., q n }, h<n, 

be a subset of the quantization levels of the hnage. Then, if u(x) G [q k , %+i], Yang’s filter is 
given by the image interpolation 

Fyu(x) = {q k+ 1 - u(x))Jq k (x) + (u(x) - q k )J<j k+1 (x), (27) 

where 

Ey e Ar(x)^(%-yK(l x -ylHy) 

E ye iv(x) £h(<ik - y)«'p(l x - y|) ' 1 j 

Thus, if h = n then Fy coincides with the exact discrete version of the bilateral filter, whereas 
if h < n an approximation of the bilateral filter is obtained by interpolation on the quantization 
range kernel levels. Observe that since x is no longer present in the range kernel of .Jq k (x), 
both factors in (]28[) may be computed by fast approximated convolution algorithms, which is 
the main strength of the algorithm provided in [39] 

However, in the first experiment (Yaroslavsky filter), and in order to obtain exact results 
when n = n, we implemented Yang’s algorithm by using our rearranged version (l24f) apphed to 
Jq k (x), explaining the increment of execution times when compared to the actual implementation 
of Yang et al. 

As already mentioned, we also use for comparison purposes the fast filtering using the per- 
mutohedral lattice introduced by Adams et al. This algorithm is based on resampling techniques 
and specially useful for high dimensional filtering. We refer to [2] for a thoroughfull explanation 
of the method. 

We use the following notation for the algorithms to be compared: 

• YPB: Yaroslavsky pixel-based version, formula. (1231) . 

• YRR: rearranged Yaroslavsky version, formula (f23|) . 

• EYangn: exact Yang’s algorithm, formula (|27[) . computed with the rearranged strategy, 
for n = 256, 64 or 8 interpolators. 

• BPB: Bilateral pixel-based version, formula (125(1 . 

• BR.RM: rearranged bilateral version, formula (f26j) . with the maximum number of spatial 
kernel discretization levels, m(4) = 42, m(8) = 135,... 

• BR.R20: rearranged bilateral version, formula (f2(Ij). with fixed number of spatial kernel 
discretization levels, m = 20. 

• Yn: Yang’s code provided by the authors, see footnote 1 in page [IS] for n = 256 or h = 8 
interpolators. 
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• P: Permutohedral Lattice bilateral filter code provided by the authors, see footnote 2 in 
page U51 

5.4 Results 

In the first experiment, we tested the approximation quality of the rearranged and the exact 
Yang’s versions of the Yaroslavsky filter to the pixel-based version. In Table [I] we show the 
PSNR of YRR and EYang with respect to YPB. We recall that, according to [32]. PSNR values 
above 40dB often correspond to ahnost invisible differences. The high PSNR. values obtained 
for YRR. and EYang'256 show that these algorithms give practically the same results than YPB, 
up to rounding errors. EYang64 still gives very good approximation results, whereas these 
significantly diminishes for EYang8. 

In Table [T] we also show the execution tunes. YRR is up to 400 times faster than YPB. 
The execution times of EYang compared to YRR. are always higher, but this has to be taken 
cautiously since the code implementations were not optimized. In this experiment we were more 
interested in approximation quality. 

In the second experiment, we tested the approximation quality of the algorithms described in 
the previous subsection to the ground truth image and to the bilateral pixel-based filter (BPB). 

In Table [2] we show the corresponding PSNR’s. We see that all the algorithms employed 
give similar results when compared to the ground truth image. Thus, if this were the choice 
criteriuin, the faster, that is Y8, should be considered. 

However, when compared to the BPB, the PSNR’s are quite different. As expected, the 
rearranged version with the maximum number of quantization levels for the spatial kernel, 
BR.R.M, gives the same result, up to rounding errors, than the BPB. However, the Yang’s 
algorithm with the maximum number of levels, Y256, although should give exact results, it 
does not, revealing other sources of error beyond rounding errors. Even the BRR20 and the 
Yaroslavsky YRR. give better approximations than Y256. In particular, BRR.20 has always 
values of PSNR. around 40dB. The Permutohedral algorithm gives similar results than Y256 and 
Y8. 

Let us mention two difficulties we found with the Yang code implementation. The first 
is related to low h— values, for which Y8 gives poor results due to some artifact formation. 
The second is related to large images, for which huge memory allocation requirements result in 
running time execution errors. 

In Table [3] we collect the execution times obtained in tins experiment. Only for the smaller 
image sizes and h —values Y8 has a competitor in YRR.. For large images and high h —values 
the second fastest algorithm is, as expected, the permutohedral, although always taking around 
twice the time of Y8. Both BRRM and BRR.20 give execution times considerably higher than 
the other algorithms, for our non-optimized codes. 

Finally, in Figures [2[ [3] and H] we show some of the filters outcome. The columns give, from 
left to right, the denoised image, a detail of the image, the contour plot of the detail, and the 
histogram of the denoised image. The rows give, from above to below, the noisy image, the 
results of BPB, BRR.20, Y8, YRR and P, and in the last row, the ground truth clean image. 

It may be observed that, in general, the permutohedral algorithm performs a hardest denot¬ 
ing than the other algorithms. Tins is also reflected in the formation of picks in the histogram. 
The Y8 gives the softest denoting and performs some smoothing in the histogram of the denoised 
image, unlike the other algorithms. Finally, BPB and BR.R.20 are ahnost indistinguishable, while 
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YRR is very close to them, although with a harder denoising effect. The staircasing effect due 
to the reduction of level sets is present in all the algorithms, although not specially relevant in 
none of them. 


6 Summary 

In this paper we introduced the use of functional rearrangements to express bilateral type filters 
in terms of integral operators in the one-dimensional space [0, |D|]. 

We have proved the equivalence between the pixel-based formulation of the original filter 
and its rearranged version. In addition, we proved the convergence of discrete finite step-wise 
approximations of image and filter to their corresponding continuous limits. Although this is 
a property easy to obtain in the pixel-based formulation, it is far from trivial in its rearranged 
version. 

In the case in which the spatial kernel, w p , is homogeneous (e.g. the Neighborhood filter), it 
was proven in [23J that, the level set structure of the image is left invariant through the filtering 
process, allowing to compute the filter jointly for all the pixels in each level set, instead of 
pixel-wise. 

However, if the spatial kernel is not homogeneous the invariance of the u —level sets through 
the filter is, in general, broken. Despite this fact, there still remains an important property of 
the rearranged version: the range kernel /C/,('u(x) — u(y)) is transformed into a pixel-independent 
kernel K-(qk. — Qi), implying a large gain in computational effort, as already observed for particular 
cases in [37] . We have illustrated numerically this property. 

The present state of the rearranging technique has an inherent restriction: the need of a 
relation of order in the range space. While this may be a minor issue for vector-valued color 
images, for which a map to some one-dimensional color space may be used, it is unclear how to 
extend the method to patch-based algorithms such as the Nonlocal Means. Whether or not the 
patch reordering techniques introduced, among others, by Ram et al. [IT] may be adapted to 
the rearranging approach will be the focus of future research. 
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h 

YPB 

YRR 

EYang256 

EYang64 

EYang8 

Clock (256 x 256) 

h 

Time 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

4 

0.28 

0.02 

61.09 

0.09 

61.09 

0.08 

37.02 

0.10 

19.47 

8 

1.02 

0.02 

59.93 

0.12 

59.93 

0.13 

42.92 

0.11 

25.01 

16 

3.91 

0.03 

55.90 

0.15 

55.90 

0.15 

44.23 

0.14 

21.84 

32 

15.40 

0.06 

49.73 

0.10 

49.73 

0.09 

43.37 

0.10 

21.36 


Boat (512 x 512) 

h 

Time 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

4 

1.10 

0.06 

66.21 

0.42 

66.21 

0.37 

36.53 

0.38 

20.64 

8 

4.23 

0.08 

61.53 

0.44 

61.53 

0.45 

41.23 

0.50 

21.60 

16 

16.79 

0.14 

57.28 

0.51 

57.28 

0.40 

42.17 

0.40 

19.90 

32 

71.15 

0.22 

49.74 

0.48 

49.74 

0.36 

43.45 

0.36 

21.82 

Airport (1024 x 1024) 

h 

Time 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

4 

4.43 

0.23 

69.65 

1.58 

69.80 

1.44 

36.43 

1.61 

19.80 

8 

17.19 

0.30 

64.02 

1.80 

64.02 

1.78 

41.58 

1.82 

22.60 

16 

70.75 

0.58 

59.14 

1.80 

59.13 

1.70 

42.72 

1.75 

20.60 

32 

291.4 

0.92 

50.62 

1.71 

50.62 

1.44 

43.23 

1.54 

19.42 

Still life (2144 x 1424) 

h 

Time 

Time 

PSNR. 

Time 

PSNR 

Time 

PSNR 

Time 

PSNR 

4 

14.07 

0.68 

64.30 

3.82 

64.29 

3.45 

38.08 

3.97 

19.58 

8 

52.71 

1.06 

59.27 

4.41 

59.27 

3.95 

43.90 

4.00 

22.24 

16 

204.0 

1.71 

54.59 

5.20 

54.59 

4.61 

45.01 

4.61 

21.28 

32 

834.8 

2.65 

48.37 

4.74 

48.37 

4.73 

44.49 

4.62 

22.61 


Table 1: Comparison between the pixel-based Yaroslavsky (YPB) filter algorithm, the rear¬ 
ranged version introduced in this article (YRR), and several instances of Yang’s exact algorithm 
(n = 256, 64, 8). For each algorithm, left column gives execution tunes (in sec.) while second 
column gives PSNR with respect to YPB. Image sizes enclosed in parentheses. 
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Compared to ground truth Compared to BPB 

Clock (256 x 256) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

3.62 

3.62 

3.62 

3.32 

-1.24 

3.58 

3.64 

66.2 

42.6 

26.9 

1.17 

25.9 

29.7 

8 

4.08 

4.08 

4.08 

3.77 

3.89 

4.34 

4.03 

62.2 

43.1 

26.5 

21.3 

19.2 

27.2 

16 

4.40 

4.40 

4.40 

4.04 

4.10 

4.30 

4.31 

63.8 

41.5 

23.2 

20.3 

13.5 

21.1 

32 

4.31 

4.31 

4.31 

3.93 

4.04 

3.31 

3.55 

57 

38.4 

16.5 

16 

7.17 

12.4 

Boat (512 x 512) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

10.60 

10.60 

10.60 

10.80 

4.36 

10.60 

10.60 

69.20 

42.70 

26.80 

5.26 

25.80 

29.40 

8 

11.00 

11.00 

11.00 

11.00 

10.50 

11.00 

10.70 

63.10 

42.10 

26.30 

20.60 

17.70 

26.00 

16 

9.23 

9.23 

9.23 

8.93 

8.80 

7.46 

8.69 

62.90 

40.70 

23.10 

18.30 

12.40 

20.30 

32 

4.83 

4.83 

4.83 

4.20 

4.22 

2.59 

3.95 

57.70 

38.40 

17.20 

15.80 

7.80 

14.20 

Airport (1024 x 1024) 

k 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BR.R20 

Y256 

Y8 

P 

YRR 

4 

3.83 

3.83 

3.83 

ND 

3.15 

3.81 

3.83 

71.00 

42.60 

ND 

11.00 

26.00 

29.40 

8 

3.89 

3.89 

3.89 

ND 

3.85 

3.78 

3.87 

67.00 

42.20 

ND 

21.70 

18.40 

26.70 

16 

3.60 

3.60 

3.59 

ND 

3.59 

2.94 

3.59 

64.30 

41.10 

ND 

17.00 

12.60 

21.70 

32 

2.45 

2.45 

2.44 

ND 

2.35 

1.07 

2.38 

58.20 

38.70 

ND 

11.5 

8.53 

15.50 

Still life (2144 x 1424) 

k 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

-1 

-1 

-1 

ND 

-1.71 

-0.99 

-0.99 

66.60 

43.30 

ND 

7.96 

25.50 

30.20 

8 

-0.89 

-0.89 

-0.89 

ND 

-1.3 

-0.83 

-0.90 

63.00 

43.60 

ND 

19.10 

19.90 

28.40 

16 

-0.92 

-0.92 

-0.92 

ND 

-1.39 

-0.97 

-0.97 

60.70 

42.80 

ND 

15.40 

16.90 

23.00 

32 

-1.31 

-1.3 

-1.3 

ND 

-1.94 

-1.47 

-1.48 

54.70 

40.50 

ND 

10.50 

12.30 

17.00 


Table 2: Left box: PSNR between the ground truth image and the pixel-based Bilateral 
filter (BPB), its rearranged version for the maximum number of meaningful levels (BRRM), the 
rearranged version for 20 levels (BRR20), two instances of Yang’s algorithm, (Y256 and Y8), the 
Permutohedral algorithm (P), and the Yaroslavsky filter in its rearranged version (YRR). Right 
box: PSNR between the pixel-based Bilateral filter (BPB) and the other algorithms. “ND” 
means “no data available “ due to too huge Y256 memory requirements for these images. Image 
sizes enclosed in parentheses. 
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Execution time Execution time relative to BPB 

Clock 

256 x 256) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

0.56 

0.23 

0.17 

0.32 

0.02 

0.46 

0.02 

2.43 

3.29 

1.75 

28 

1.217 

28 

8 

2.06 

0.52 

0.25 

0.3 

0.02 

0.19 

0.02 

3.96 

8.24 

6.87 

103 

10.84 

103 

16 

7.98 

2.56 

0.62 

0.3 

0.02 

0.05 

0.03 

3.12 

12.9 

26.6 

399 

159.6 

266 

32 

31.41 

11.4 

1.8 

0.32 

0.02 

0.03 

0.06 

2.76 

17.4 

98.2 

1570 

1047 

524 

Boat (512 x 512) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

2.19 

0.9 

0.68 

1.36 

0.1 

2.11 

0.06 

2.43 

3.22 

1.61 

21.9 

1.038 

36.5 

8 

8.33 

2.17 

1.12 

1.24 

0.06 

1.04 

0.08 

3.84 

7.44 

6.72 

138.8 

8.01 

104 

16 

32.64 

10.1 

2.44 

1.2 

0.06 

0.32 

0.14 

3.23 

13.4 

27.2 

544 

102 

233 

32 

132.1 

47.2 

7.27 

1.16 

0.06 

0.14 

0.22 

2.79 

18.2 

114 

2201 

943.2 

600 

Airport 

1024 x 1024) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

8.76 

3.56 

2.78 

ND 

0.24 

9.88 

0.23 

2.46 

3.15 

ND 

36.5 

0.8866 

38.1 

8 

33.37 

8.45 

4.29 

ND 

0.22 

5.08 

0.3 

3.95 

7.78 

ND 

151.7 

6.569 

111 

16 

132.8 

39.6 

9.36 

ND 

0.22 

1.37 

0.58 

3.35 

14.2 

ND 

603.8 

96.96 

229 

32 

529.5 

185 

28.1 

ND 

0.2 

0.58 

0.92 

2.86 

18.8 

ND 

2648 

912.9 

576 

Still life 

2144 x 1424) 

h 

BPB 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

BRRM 

BRR20 

Y256 

Y8 

P 

YRR 

4 

26.5 

10.3 

8.02 

ND 

0.64 

31.4 

0.68 

2.58 

3.3 

ND 

41.41 

0.8429 

39 

8 

99.07 

24.7 

12.2 

ND 

0.62 

9.35 

1.06 

4.01 

8.09 

ND 

159.8 

10.6 

93.5 

16 

384.8 

111 

26.8 

ND 

0.66 

2.58 

1.71 

3.47 

14.4 

ND 

583 

149.1 

225 

32 

1536 

511 

80.6 

ND 

0.62 

1.52 

2.65 

3 

19.1 

ND 

2478 

1011 

580 


Table 3: Execution times for the images used in the experiments (sizes enclosed in parentheses). 
Left box: Execution times ot the pixel-based Bilateral filter (BPB), its rearranged version for the 
maximum number of meaningful levels (BRRM), the rearranged version for 20 levels (BRR20), 
two instances of Yang’s algorithm, (Y256 and Y8), the Permutohedral algorithm (P), and the 
Yaroslavsky filter fix its rearranged version (YRR). Right box: Execution times relative to BPB, 
that is, number of times the other algorithms are faster than the BPB algorithm. “ND” means 
“no data available “ due to too huge Y256 memory requirements for these images. 
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Figure 2: Image Clock, h = 8. The columns are, from left to right, the denoised image, a detail 
of the image, the contour plot of the detail, and the histogram of the denoised image. The rows 
give, from above to below, the noisy image, tlfi?results of BPB, BRR20, Y8, YRR and P, and 
in the last row, the ground truth clean image. 
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Figure 3: Image Boat , h = 16. The columns are, from left to right, the denoised image, a detail 
of the image, the contour plot of the detail, and the histogram of the denoised image. The rows 
give, from above to below, the noisy image, tH&3results of BPB, BRR20, Y8, YRR and P, and 
in the last row, the ground truth clean image. 
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Figure 4: Image Still life, h = 32. The columns are, from left to right, the denoised image, a 
detail of the image, the contour plot of the detail, and the histogram of the denoised image. The 
rows give, from above to below, the noisy ima$fl the results of BPB, BRR20, Y8, YRR and P, 
and in the last row, the ground truth clean image. 

































































